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Various phenomenological models of particle multiplicity distributions are discussed using a gen- 
eral form of a unified model which is based on the grand canonical partition function and Feynman's 
path integral approach to statistical processes. These models can be written as special cases of a 
more general distribution which has three control parameters which are a, x, z. The relation to 
these parameters to various physical quantities are discussed. A connection of the parameter a with 
Fisher's critical exponent r is developed. Using this grand canonical approach, moments, cumu- 
lants and combinants are discussed and a physical interpretation of the combinants are given and 
their behavior connected to the critical exponent t. Various physical phenomena such as hierarchical 
structure, void scaling relations, KNO scaling features, clan variables, and branching laws are shown 
^ in terms of this general approach. Several of these features which were previously developed in terms 

, of the negative binomial distribution are found to be more general. Both hierarchical structure and 

' void scaling relations depend on the Fisher exponent r. Applications of our approach to the charged 

particle multiplicity distribution in jets of L3 and HI data are given. It is shown that just looking at 
the mean and fluctuation of data is not enough to distinguish these distributions or the underlying 
mechanism. The mean, fluctuation and third cummulant of distribution determine three parameters 
^ , X, z, a. We find that a generalized random work model fits the data better than the widely used 

' negative binomial model. 

o 

o 

m 

^ I. INTRODUCTION 

Pion multiplicity distribution and their associated fluctuation and correlations have been of interest for several 
O reasons. Several models predict large fluctuations such as the disoriented chiral condensate [1,2] and in density 
2 fluctuations from droplets arising in a first order phase transition [3] . A well known procedure for studying correlations 
uses the Bose-Einstein symmetries associated with pions in a Hanbury Brown- Twiss analysis [4]. Such an analysis 
gives information about the space time history of the collision through measurements of source parameters. If the 
density of pions becomes large, Bose-Einstein correlation may also lead to a strongly emitting system which has been 
called a pion laser [5]. The pion laser model has been recently solved analytically by T. Csorgo and J. Zimanyi [6]. 
' The importance of Bose-Einstein correlations has also been illustrated in the observation of a condensation of atoms 
in a harmonic oscillator or laser trap [7] . Previous interest in pionic distributions have centered around the possibility 
of intermittency behavior [8] and fractal structure based on parallels with turbulent flow in fluids. A distribution 
widely used to discuss such features has been the negative binomial (NB) distribution [9] with its associated clan 
structure [10,11] and KNO scaling feature [12]. KNO scaling properties have been interpreted in terms of a phase 
transition associated with a Fcynman- Wilson gas [13]. Various other issues associated with pions include evidence 
for thermalization [14], critical point fluctuations [15,16], fluctuations from a first order phase transition [17], charge 
particle ratios and question of chemical equilibrium [18], the behavior of fiuctuations in net charge in a QG plasma 
for transition [19,20]. 

For lower energy heavy ion collisions, multifragmentation of nuclei takes place. The fragment distribution can also 
be described statistically by considering all the possible partition of A nucleons into smaller clusters [21,22]. This study 
gives a tool for the description of nuclear multifragmentation distributions [23] , nuclear liquid-gas phase transition 
[24], critical exponent, intermittency, and chaotic behavior [25,26] of nuclear multifragmentation. The same model 
can describe pionic distribution. This possibility arises in our approach which is based on Feynman path integral 
methods where symmetrization of A bosons or anti-symmetrization of A fermions leads to a cycle class decomposition 
of the permutations associated with these symmetries. The correspondance comes from the identification of clusters 
of size k and the cycles of length A: in a permutation as discussed below. In nuclear multifragmentation the number 
of nucleons. A, is fixed and thus a canonical partition function approach is appropriate. By contrast, the number of 
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pions, A, in particle production is not fixed and thus wc need to use a grand canonical partition. Using this parallel, 
we also note that some results from cluster yields can be carried over into particle multiplicity distributions and 
associated properties. Namely, we show the importance of the Fisher critical exponent r and relate it to one of the 
parameters called a in our approach which has three main parameters a, x, z. Moreover, taking special values of a, 
or equivalently r, reduces our unified model to various specific cases that are frequently used in particle production 
phenomenology. Quantities that appear in this development can also be related conibinants as will be discussed. In 
so doing we can give a physical significance to the combinants and show how the Fisher exponent appears in them 
and how the resulting hierarchical structure and void scaling relations also depend on its value. 

In next section a summary of a generalized statistical model based on a grand canonical partition function will 
be given. This generalized model can be used to study multiplicity distributions associated with particle production 
such as at RHIC. Various models of particle multiplicity distribution using the grand canonical partition function 
will then be developed. In Section II-B, we will summarize the multiplicity distribution and its various moments 
of all order for a general grand canonical ensemble which we will then use for particle production. The physical 
meaning of the parameters used in this model is given here. The differences and relations between the canonical 
nuclear multifragmentation and the grand canonical multiparticle production are sumarized in Sect. II-C. In Sect. 
II-D and E, moments (especially the mean and the fluctuation) in the general grand canonical model are related to 
variables used in other existing standard models in describing multiparticle production. The connection of the general 
grand canonical model developed here to existing standard models are summarized in Sect. II-F. Moreover, we derive 
a generalized model (HGa) of the grand canonical partition which can further be reduced to a geometric, negative 
binomial, and Lorentz/Catalan model. Various moments and multiplicity distribution for a generalized HGa model 
are summarized in Sect. II-G and another generalized model (GRWID) in Sect. II-H. In Sect. Ill, we compare various 
statistical properties between different models within the generalized HGa model and with data for charged particle 
distribution in jets of L3 and HI. Fluctuation and multiplicity distribution in HGa are compared with other standard 
models of multiparticle production. 

II. GENERALIZED PROBABILITY DISTRIBUTIONS 

Consider a system composed of N different types of species or objects which could be the fragments in a 
fragmentation or in a cycle class description. Any event of such a system can be associated with a vector 
n = {uk} = {ni,n2, ■ ■ ■ ,nk, ■ ■ ■ jUn) or l"i2"^3"'* • • • fc"*^ • • • A^"" where the non-negative integer Uk is the num- 
ber of individuals of species k. For example can be the number of clusters of size k or the number of cycles of 
length in a given permutation of n particles. The later is important for Bose-Einstein and Fermi-Dirac statistics 
and particle multiplicity distributions. A general block picture of n is shown in Fig. la. Fig. lb shows how the various 
partition can be developed as an evolution from successively smaller systems. The number of species A'^ can be infinity 
in general. 




FIG. 1. Building partitions with blocks. 
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Various probability distributions related with this system can be developed by assigning an appropriate weight Xk 
to each type k. A weight W{x, n) is then given to each event n and the type of weight that will be considered has 
the structure: 



N r 



fe=i 



(1) 



The rife! are Gibbs factorials. Such a weight structure appears in Fcynman's path integral approach and Bose-Einstein 
problems [27]. The X]^ will be given below and contains various physical quantities. In Sect. II B wc will show that the 
XkS are also the combinants which in turn can be related to the factorial cumulants. Summing the weight W{x,n) 
over all the possible events of n, the grand canonical partition function Z{x) of the system is given as 



N 



Z{x)=Y,W{x,n) = Y,\[ 



k=l 



exp 



■ N 



Lk=l 



(2) 



The last equation holds due to the form of Eq.(l) of the weight W{x,n), i.e., the factor x^'=/nfc! is the ri/j-th order 
expansion term of e^'' . 

Introducing other quantities ak to each individual entity or group of type k, a = {cKfe} = (en, 0L2, ■ ■ ■ , cx-ki " " " ) oln), 
we can define a canonical partition function Za{x) with a fixed A as 



JV 



A = ^afe/ifc = a-n 
fe=i 



Za{S) = Y,W{ 



N 



x, nj 



En 

UA k=l 



(3) 
(4) 



with Z{x) = Za{x). Here X^^^ is the summation over all events with a fixed value of A, i.e., over a canonical 
ensemble of a fixed A, and the J2a ^ summation over all the possible values of A; it becomes Y1a=o '^^^ °^ 

ak = k with positive integer k. The multiplicity of a partitions is 
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M = J2nk 



(5) 
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The case = A; is encountered in fragmentation problems and permutation problems. Various type of can be 
used depending on the physics related with the quantity A as discussed in Ref. [22]. The physics of the canonical 
ensemble depends on the choice of the quantity ak [21,22]. If we take Xk cx z"'', then the canonical partition function 
Za{x) is the dependent term of the grand canonical partition function Z{x). If we take Xk oc x, then x counts 
the number of clusters M explicitly and the dependent term of the grand canonical partition function Z{x) is 
the partition function for events with multiplicity M . There always is at least one event having A = 0, i.e., the event 
where all n/c's are zero, n = 0. Thus Zq{x) = 1 if all are non zero positive since then there are no other possible 
events having A = 0. Due to the form of the weight W{x, n) given by Eq.(l) the canonical partition function Za{x) 
satisfies a recurrence relation [25,27] 



Za{x) = ^ E '^kXkZA- 



{x) 



(6) 



This relation is nothing but the constraint Eq.(3) in terms of the mean {nk)A using Eq.(8). For non-zero positive ak, 
there is no case having A < 0, i.e., = for A < and thus the Za can be obtained by the recurence relation of 
Eq.(6) starting from Zo{x) = 1. 



A. Probability distribution 



In a canonical ensemble of fixed A, we can define a probability distribution of a specific partition n as 

W{x, n) 



PA{x,n) 



Za{x) 



(7) 
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With this probabihty, various mean vahies, fluctuations and correlations of the number of species rifc can be evaluated 
as a ratio of canonical partition functions for two different values of A such as [21,22] 



(rife - m)! {iij - ly. 



4^ (rife - m)\ {rij - l)\ 



PA{x,n) =x^x\ 



Za- 



■maf^ — L 



Za{x) 



(8) 



Thus we have {nk)A — x^ ■ The recurrence relation of Eq.(6) then follows simply using the fact that 

A = Ylik=i '^k{nk)- This distribution has been used in describing various fragment distributions in nuclear multi- 
fragmentation [21-23]. 

Now knowing the partition functions Z{x) and Za{x), we can associate a probability Pa{x) of the system to have 
a fixed value of A in a grand canonical ensemble as 



Pa{x) 



Za{x) ^ 1 1 
Z{x) Z{x) r{A + 1) 



Z{x, z) 



2=0 



Z{x, z) = Za{x)z^ = exp 



E 



XkZ 



(9) 
(10) 



The last step follows from the fact that the z"^ power term of Z{x,z) is Za{x) if we put Xk = z°"'. For particle 
production Pa{x) is the probability of having A particles. Thus the generating function Z{x,z) of Pa can also be 
looked at as a grand canonical partition function with the weight Xk replaced to be Xkz"'' ^ where the variable z counts 
A explicitly and Z{x, 1) = Z{x). If we consider Pa{x, z) = Pa{x)z^ then z has two roles; one as a weight which is 
assigned the same to each constituent and another as a generating parameter of the probability Pa{x). For the case 
that Zq(x) = 1, the void probability Pq [28] is the inverse of the grand canonical partition function, i.e.. 



(11) 



Inversely the canonical partition function Za{x) is the probability Pa normalized or rescaled by the void probability 
Po, i.e., 



Za{x) = 



r(^ + i) 



Z{x, z) 



Pa{S) 
Po{S) 



(12) 



Another type of generating function of Pa that is frequently used may be defined as 

- oo 

G{x,u) = j2PA{m - u)^ = ^ E^^(^)(i - 



Z{x, z = 1 — u) 



Z{x) 



exp 



Y^xkKi-ur"-!] 



Pa{x) 



1 



r(^ + i) 



--) G{x,u) 



(13) 
(14) 



We see that G(f , 0) = 1, G{x, 1) = Po{x) = Zo{x)/Z{x) and G{x, 1- z) = Z{x, z)/Z{x). Once the probability Pa{x) 
is determined, various statistical quantities can be evaluated. Also the grand canonical partition function Z{x, z) is 
given once we know the thermodynamic grand potential 



fi(x, z) = — \nZ{x, z) = — y~^a:fc^° 



(15) 



This result can be used to study the statistical properties of a system. Moreover various moments and cumulants, 
mean values and fluctuations may be obtained using the generating function [21,22]. 
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B. Moments and cumulants; combinants and hierachical structure in grand canonical ensemble 



This subsection gives general expression for various quantities that will be used later when we discuss specific 
models. Since the probability of a specific event n in a grand canonical ensemble is given by 



(16) 



the mean of a quantity P in a grand canonical ensemble is related to the mean of F in a canonical ensemble as 

{F) = J2 FP% n)=J2 P^i^) E ™(^' ^) = E Pa{x){F)a (17) 



We can easily show that 



('^fe) = E "^kPix, «) = E PA{x){nk)A = Xk 

n A 

(M) = ^ ( ^ nfc J P{x, n) = Y,Xk = -0(f , z = I) 

fi \ k / k 



(18) 
(19) 



This result shows that the weight factor Xk in this model is the mean number (n/t) in a grand canonical ensemble 
and that the mean multiplicity (M) is related to the thermodynamic potential. The m-th power moment of A and 
its factorial moments are given simply by 



n \fc=l 



r(^ + i) 



r{A + i) 



-P{x, n) 



Z(x) 
1 



Z{x,z) 



-) Zix,z) 



z = l 



^ T{A - m + 1) / ' ' ^ r(A - TO +1) ' ' ' Z{x) 
Similarly the m-th cumulants , which is the power moments of , and the factoral cumulants are 

d 



(a™)(f) 

fm{x) 



E""«^) = 

\k=l I 



dz 



In Z{x, z) 



= E "^kXk 



(20) 
(21) 

(22) 



E 



r(afe-TO + l) 



rik 



-) lnZ(f,.) 



d_ 
du 



\nG{x,u) 



N 



r(afc + l) 
„=o ^r(afe-m + l)''^'= 



Due to Eqs.(lO) and (15) we can see easily that, for the power moments of ak, 



N 



(aO)=^Xfc=/o = (M) 
fe=i 

TV 

("^) = E = /i = (^) 
fe=i 

N 

J2 ^l^^ =f2 + fi = {{A {A)r) = {A') (A)' = a' 



fe=i 

N 



{a') = J2 <^kXk =h + 3/2 + /i = {{A - {A)f) 



(23) 

(24) 
(25) 
(26) 
(27) 



fe=i 



The power moments of are directly related to the power moments of A measured from the mean (A), i.e., the 
cumulants (a™) are same with the central moments of A. This simple relation does not hold for to > 4 but we can 
evaluate them starting from (a^) using the recurrence relation 
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, \ m+l / , 



(28) 



Similarly the m-th factorial cumulants fm, which is the factorial moments of ak, can be found using recurrence relation 



fm+l{x,z) 



z ( ^ ) ln/m(x, z) - m or 



fm 

{x,z) _ d ( fjn-i{x,z) 



dz 



starting from 



yvn—l 



N 



Tz 



fe=i 



Mx,z) = {M){x,z) = \nZ{x,z) = -fl{x,z) = ^x^^"'' 

N 

fi{x, z) = {A){x, z)=^ akXkz""- 



(29) 

(30) 
(31) 



fe=i 



The reduced factorial cumulants Km defined in Ref. [28] corresponds to the factorial cumulants /„ normalized with 
mean number {A) = A as 



^m{Xi Z^ 



fm{X-, Z^ 



(32) 



Thus with Ki = 1 and kq = fo = {M) . The cummulants are directly related with the generating functions as 

G{x,u)=exp Vl^/„(f) =exp V Mi_!^(«™)(f) (33) 

_m=l J \_m=l 

Zix,z) = e^p V =exp V iiif^(a™)(f) (34) 



.m=0 



m! 



.m=0 



m! 



The generating function Z{x, z) differs from the generating function G{x, u = 1 — z) only by an extra term of m = 

in their exponent. 

The relation between the x^'s and the Z shows that the x^'s arc also the conibinants of Ref. [29]. In the approach 
presented here, the combinants are given an underlying significance through the partition weight W{x,n) of Eq.(l). 
In turn the combinants Xk can be related to the factorial cumulants fm defined by 



In Z{x, z) = Y^ Xkz""- = 



m! 



■fmix) 



(35) 



k m=0 

The factorial cumulants fm are the m-th order factorial moments of ak of Eq.(23). Thus 



/™ = m! ^ ( J Xfe 



(36) 



for ak = k which is the Eq.(23) with TV 00. The normalized factorial cumulant, i.e., the reduced cumulant, is 

Km = fm/A^ = {m-1)\K^-' (37) 

for a negative binomial (NB) distribution, This result of Eq.(37) shows that for NB has an hierarchical structure 
of a distribution at the reduced cumulant level which was realized for the NB distribution in Ref. [28]. This result 
will be generalized later. 

Also using the above power moments and factorial moments we can study voids and void scaling relation, hierarchical 
structure, combinant and cummulant properties which will be discussed below. 
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C. Multi-fragmentation versus multiparticle production and the Fisher exponent for pcirticle production 

Since our approach was first used to discuss multifragmentation and then later extended to include multiparticle 

production, wc briefly mention some of difference between multifragmentation and multiparticle production. We also 
show how the Fisher exponent r, which initially appeared in cluster yields around a crititical point, may manifest 
itself in particle production yields. 

In nuclear multifragmentation, ak = k is the number of nucleons in a fragment and rife is the number of fragments 
of size k. The total number of fragments is M = and the total number of nucleons is ^ = X)fe ^'^k- In nuclear 

multifragmentation, the total number of nucleons A is usually fixed and we study the distribution in size k of the mean 
multiplicity (nfc)^ = XkZA-k{x)/ZA{x) of fragments of size fc in a canonical ensemble. In multiparticle production the 
A is the total number (multiplicity) of produced particles and is not fixed and we study the multiplicity distribution 
of produced particles Pa = Za{x)/Z{x) as a function of A in a grand canonical ensemble. The mean multiplicity 
{M)a and {nk)A in a canonical ensemble are related with grand canonical ensemble as 

oo 

{rik) = ^^{nk)APA = Xk 

oo 

{M) = J2{M)APA = Y.Xk 

A=0 k 
oo oo 

= $Z 5Z Hnk)APA = ^ ^Pa = 5Z ^^fe 

A=0 k A=0 k 

The weight Xk can be determined experimentally from the mean multiplicity (rifc) of cluster size /c in a grand canonical 
ensemble. 

In our initial fragmentation studies [21-23] we have used Xk = xz^ jk. This choice gives Za{x,z) = (z^/A!)r(x + 
A)lV{x) and thus ZA-kjZA z~'' as ^ ^ oo where the z = e^^ with fj, the chemical potential. From Eq.(8), 
{nk)A = XkZA^k/ZA, thus {nk)A — > x/k as A — > oo in the case of nuclear fragmentation in the canonical ensemble 
of fixed A for Xk = xz^ /k and the yields fall as a power law at a critical point (x = 1) for this choice of Xk- 
The yield in the grand canonical ensemble is (n/c) = Xk = xz^ jk and this becomes the same as ijik) a x/k in 
A — > oo limit for ^ = 1 with /i = 0. In general, cluster yields fall for large A; as a power law l/fc"^ at a critical point 
where r is the Fisher critical exponent. The above choice of Xk can be generalized to Xk = xz^ /k'^ so that the cluster 
yields fall as 1/k'^. The r determines the grand canonical partition fimction Z{x) = exp[^j,Xfc] = exp[x^j, 2;'^/fc'^]. 
Forms for Xk whose asymptotic behavior is xz^ /k'^ can also be used. In section II. F we will give various choice 
for Xk for particle production. For example, the above choice of Xk with r = 1 gives the frequently used negative 
binomial distribution with x the negative binomial parameter which determines the degree of depature from Poisson 
statistics. Other distributions can be generated and their associated r dependence will also be given. Then we will 
show how the r dependence of the combinants shows up in hierarchical structure and void scaling relations in sect 
II. G and in other places. The physical meaning of x and z will also be discussed for these other cases. Since an 
exact description of particle multiplicity yields is not known, we consider t as a free exponent to be determined 
by comparing various resulting distributions with experiment. This procedure is what is done in cluster yields. In 
multi-fragmentation distribution, the mean number of clusters (M) is determined by x and z is directly related with 
the number of constituent particles A. The large k behavior of the mean number of cluster (rife) gives the Fisher 
power T. In multiparticle distribution, x and z determine the mean particle number {A) or the peak position of the 
particle number distribution Pa and the fluctuation = (A^) — (A)"^ or the width of the particle distribution. The 
Fisher power t determines the large k behavior of the weight Xk and offers a new parameter in particle production, 
(nfc) = Xk in grand canonical ensemble. 

Bose-Einstein condensation of atoms in a box of sides L of dimensions d have x = L'^/Ay with Xt ~ h/{2TrmkBTy^'^ 
and have r = 1 + d/2. Feynman used random walk arguments, the closing of a cycle parallels a closed random walk, 
to discuss his choice of Xk in his discussion of a superfluid phase transition in liquid helium. Thermal emission 
of pions based on statistical mechanics and equilibrium ideas have been popular descriptions of pions coming from 
relativistic heavy ion collisions. For thermal models [30], the Xk = {VT^ /2'!T'^)(rn/T)'^ K^ikm/T) / k for a cycle length 
fc or a cluster of size k with K2 a Mac Donald function. For low temperatures, Xk = {y/)^){e~™'/'^)^/k^/'^ and the 
Boltzmann factor in mass, e"*^™/-^, suppresses large fluctuations. In the high temperature limit and/or zero pion 
mass limit Xk = {V/7r'^)T^ /k^. Thus the high temperature limit has r = 4 and the low temperature limit has t = 5/2. 
The Xk can be used to generate the pion probability distribution P„. The thermal models can be combined with 
hydrodynamic descriptions and an application was given [30] to 158A GeV Pb+Pb data measured by the CERN/NA44 
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and CERN/NA49 collaborations. The results of Ref. [30] showed a Gaussian distribution with a width about 20 % 
larger than the Poisson result. 



D. Clan parameters and void parameters and void scaling relations 



Van Hove and Giovannini have introduced clan variables Nc and ric to describe a general class of probability 
distributions, with most discussions of these variables centering around the negative binomial distribution [11]. These 
variables are defined as 

N, = {M)=\nZ = fo, rie = {A)/N, = f,/fo (38) 

where the mean number of clans is Nc and the ric is the mean number of members per clan. The Z is the grand canonical 

generating function and thus Nc = X]fe ^'k where X). is the cycle class weight distribution x. The (A) = k{nk) is 
the mean number of total members (particles). The n/j here is the number of clans of size k having k members. 

The clan variable Nc is also related to the void probability Pq = Zq/Z = 1/Z = e~^'; thus A^c = — InPo = fo- An 
important function in void analysis is x = — InPo/ (^4) = Nc/ {A) = Yjuc- Thus the void parameters, void probability 
Po and void function v = x [28], are equivalent to the generalized clan parameters Nc and Uc with the equivalence 
given by: 

Ux) = In Z{x) = - In Po{x) = Nc (39) 

X{x)^^=^ = nc-' (40) 

The —X is the normalized grand potential = — /o for the mean (A). 

Void analysis looks for scaling properties associated with specifically, x is a function of the combination (A)^ 
where ^ is the coefficient of {A)^ in the fluctuation cr^ = {A) + ^{A)^. Since /2 = (a^) - {a) = {{A - {A)f) - {A), 
the variance of A in a grand canonical ensemble becomes 

a^^{A')-{Af = {a')=J2alx, 

k 

= {A) + h = {A) + i{A)^ = {A) [1 + m{A)] (41) 

with ^ = K,2i i.e., the normalized factorial cumulant. Since the variance for a Poissonian distribution is the same as 
the mean, cr^ = (A), the ^ = 0; thus, the parameter £^{A) represents a degree of departure from Poissonian fluctuation 
normalized by mean {A) of the distribution. A well known non-Poissonian example is a NB distribution which has 
^ = ^ and this becomes Plank distribution with x — 1. Thus the variable x in Xk determines the strength of the non- 
Poissonian fluctuation term for the negative binomial distribution. In general, x and z determine the mean particle 
number or the peak position of the particle multiplicity distribution and the fluctuation or the width of the particle 
distribution. Using the recurrence relation Eq.(29) for /„, we can show that 

^{x,z){A){x) = x-\x,z) - z {j-^ \nx{x,z) - 1 = ^""^^1^'^^ - 1 (42) 

A NB distribution has x — lii(l + '?(^))/('?(^)) while the Lorcntz/Catalan (LC) distribution discussed in Ref. [31] 
and below has x = (■\/2^(74) + 1 — l)/(^(74)). We will study in Sect. Ill A, x ^(^)) i-®-, the void or clan variable 
vs the fluctuation for various choices of X}^ summarized in Table I. 



E. Ancestral or evolutionary variables 



The LC model, which has 

a;, = 12-2(^-1) (2(fc-^l)^^^fe^^^^ 

shown to be a useful model for discussing an under- 
lying splitting dynamics when ancestral or evolutionary variables p and /3 were introduced into x and z as discussed 
in Ref. [31]. Specifically x = 13/Ap and z = 4p(l - p) giving Xk = l3Ckp'^^~^\l — p)^ where Ck = \ {^^^_ ^ • 
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Percolation or splitting dynamics with a branching probability p and survival probability (1 — p) has a hierarchical 
topology as shown in Fig. 2. Weighting each diagram by Xk = pCkP^^^{l — p)'^, the evolutionary or ancestral variables 
are related to the clan variables Nc = (M) and ric = {A)/Nc. Ck is the munbcr of diagrams of size k shown in Fig. 2. 
For this case the evolutionary dynamics is just that of the LC model. Then with (3 set equal to 1, Xi = (1 — p), 

-p)^, etc. The interpretation of this set of x^s reads as follows: xi has 
and one diagram (Ci = 1). X2 has 1 branch (p^) leading to 2 surviving 



X2 = p(l - pf, X3 = 2p2(l - pf^ X4 = 5p3(l 
1 surviving line without a branch — p)^ 
lines ((1 — p)"^) and one diagram (C2 = 1) 



.T3 has 2 branch points 

„3- 



(jr 



leading 3 surviving lines ((1 

^4^ 



p)^) and two 

diagrams (C3 = 2), X4 has 3 branch points(]3''), 4 surviving lines ((1 — p)^) and 5 diagrams (C4 = 5), etc. In these 
evolutionary/ancestral variables the /o = Yl^k = 2x(l — ^/l — z) which determines Z is fo = /3 for all p < 1/2. 
For p > 1/2, /o = ^fe i'' longer a constant and is /o = /3(1 — p)/p- To keep /o = X^fc-^fc ^ constant without 
changing /i — (n) = J2k^(^k) = P \i-2p\ ' ^ "^oc ~ "^^^o introduced in Ref. [31]. The (/)oo = for p < 1/2 and 
(poo = /3{2p — l)/p for p> 1/2. For p> 1/2 there is a finite probability ((/)oo 7^ 0) that the splitting will go on forever 
(fc ^ cxd). In percolation above a certain p an infinite cluster is formed and ^00 is similar to the strength of the 
infinite cluster. Moreover the sudden appearance of (j)oo is similar to the behavior of an order parameter in a phase 
transition. The appearance of (j)oo can also be interpreted as the sudden appearance of a jet without pion (fc = 0), if 
we introduce xo = ^00 instead of Xoo = (poo ■ 



Since the clan variables are Nc = (M) 
or ancestral varables f3 and p, 



fo and He = {A)/Nc = /1//0, then, for the LC model with evolutionary 



(M) = fo = P 



{A)/N, 



l-|l-2p| 
2p 
2p{l-p) 



{A) =13 



^=2 



|l-2p| 

It 



|l-2p|(l-|l-2p|) 



1 



2nc-l 
lie - (1 ± l)/2 
ric - 1 



These can be reduced to = (3,71^ = (1 -p)/{l - 2p), and 2p=l- l/(2nc - 1) for p < 1/2 while = P{l-p)/p, 
ric = p/(2p— 1), and 2p = 1 + l/(2nc — 1) for p> 1/2. Since the branching probability p varies in < p < 1, the clan 
variable ric has ric > 1 with ric = 1 at p = and 1, only one member per clan in average, and ric = 00 at p = 1/2, 
infinitly many members per clan. The LC model thus connects the clan variable Uc to the probability p of branching 
in the evolutionary or ancestral picture of Fig. 2 or in a percolation model. For Poisson processes p = (no branching 
), Xk = (36ki (only unit cycles and no BE correlations) and rig = 1 (one member in each clan in average). 



1-p p(l-p)^ 
unit cycles period 2 cycles 



P 



P 



2p2(l-p)3 

period 3 cycles 




P 




5p^{l-p)'^ 
period 4 cycles 



FIG. 2. Evolutionary lines of descent in a hierarchical topology. Each branch increases the cycle length with probability p, 
survival 1 — p. The probability distribution evolves from Poisson to chaotic. For clusters each branch generates a bigger cluster. 



A connection of the LC model can also be made with a Ginzburg-Landau theory of phase transitions and a 
Feynman- Wilson Gas. These connections are discussed in Ref. [32]. where x and z are related to coefficients in the 
Ginzburg-Landau approach. 
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F. Unifying various distributions with Gauss hypergeometric series 



Once we identify an appropriate Xk for a physical system, then we may use our general model to study the statistical 

behavior of the system. The various models used in pion distribution can be related to our general model with ak = k 
by choosing Xk as a term in a Gauss hypergeometric series F[a, b; c; z); 

F{a,b;c;z)=f:^^^^'-^ (44) 
^ c,„ ml 

m— 

. , r(a + m) T(a + n) r(a + n + m — n) , ^ , ^ 

['^l" = -^ioT = r(a) na + n) = + ^''^ 

Usefull values of [a]„ are [l/2]„ = (2n)!/(n!2^"), [1]„ = n\, [2]„ = (n + 1)!. Considering only positive k, we choose 

[a]k-i[b]k^i z'' r(a + fc-l) r(c) r(6 + fc-l) 

'^'^ [c]k-i ^ r(a) r(c + fc-l) Tib) {k-l)\ ^ ' 

For this case, the thermodynamic grand potential or the generating function is 

oo 

fo{x) = fo{x, z) = log Z{x, z) = -Q{x, z) = ^Xk = xzF{a, b; c; z) (47) 

fe=i 

If we allow jets without a pion, then we may allow fc = also. For such a case, 

^ _ [aUb]k z'' _ r{a + k) r(c) T{b + k) z'^ 

[c]k k\~ T{a) T{c+k) T{b) k\ ^ ' 

oo 

fo{x) = fo{x, z) = log Z{x, z) = -Cl{x, z) = '^Xk = xF{a, b; c; z) (49) 

fe=o 

Wc can sec the only difference of the generating functions between the above two cases is the extra factor z for the 
grand potential. We will mostly concentrate on the first case, i.e., k ^ 0. 
Using Eq. (23) or the recurrence relation Eq.(29) and Eqs.(44) - (47), 

-J n{x,z)=z"'( — \ \ogZix,z) = z^( — \ xzF{a,b;c;z) 

= xI^Mziiz-+iF(a + m,b + m;c + m;z) 

[C\m 

+a;mi^%^ii^i^^'"F(a + m - l,b + m - l;c + m - 1; z) (50) 

[C\m-1 

The second order normalized factorial cumulant £, = K2 and the void variable x are then 

f2{x,z) f2{x,z) 



ax,z){A) 



{A) h{x,z) 
^-^^z^F{a + 2,b + 2;c + 2;z) + 2S^zF{a + 1, 6 + 1; c + 1; ^) 



fzF{a + 1, 6 + 1; c + 1; 0) + F{a, b; c; z) 

fo{x,z) _ .fo{x,z) 
(A) ~ h{x,z) 

F{a, b; c; z) 



fzF{a + 1, 6 + 1; c + 1; ^;) + F{a, b; c; z) 
For some values of a, b, and c, the hypergeometric function become a simple function; 



(51) 



(52) 
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F{a,b;b;z) = {1- z)-" 
F{a,b;a;z) = (1 - z)''' 

1 - (1 - zy-" 



F{a,l;2;z) 



z{l - a) 



F{1, 1; 2; z) = lim F{a, 1; 2; z) = -^-^ — ^ 

a^l Z 



(53) 



Various models of pion distributions can be related with these functions as listed in Table I. 



TABLE I. Various models with specific ciioice of Qfe = fc and Xk in liypergeometric scries F{a, b; c; z) of Eq.(44). Here fc = 
is not inc:hidcd. and tlins /o = lu Z = ^ Xk = xzF{a. h: c: z). Uvrv 1 < k < N witli N OC' (>x(;opt for Poisson wliicii has a 

finite A'.i;. fisher exponent r for each ,r;, as discussed in Seetdf ( ' are given too. 



Model 



foix)=lnZ 



Poisson (P) 
Geometric (Geo) 
Negative Binomial (NB) 
Signal/Noise (SN) 

Lorentz/Catalan (LC) 

Hypergeometric (HGa) 

Random Walk-Id (RWID) 

Random Walk-2d (RW2D) 

Generalized RWID (GRWID) 
Generalized RW2D (GRW2D) 



Nx5k,i 01 X tor k = 1, 2, ■ ■ ■ , N 

k 



-2(fe-i) (^2(fc-^f) l^^fe 

,.-.fc 

-2(fc-i) / 2(A; - i) \ 
I k-1 



-2{fc-l) / 2(fc- i) 
[°h-l 



iVa; = 

1-z 

—a; ln(i — 2) 

^^-x\n{l-z) 

2x\l - (i- z)'/'] 

xzF(\,\-\;z) 
xzil-z)-" 
xzF{a,a; l;z) 



1 
1 

1 
i 

1(6) 

1 
2 

b 

a 



2 
2 
1(6) 

1 
6 
1 





1 

3/2 
2 - a 
1/2 

1 

1-0 
2(1 - a) 



TABLE II. Factorial cumulants for vaxious choices of xt of Table. I 



/o =logZ 



/i = (A) 



Z-^frr 



Nx = A 

Xj^ 

1 — z 

— xln(l — z) 
^^-x\n{l-z) 
2x[l - (I - zf'''\ 



-(l-z) 



T72 



a; 



1 — 2 

X- 



a; 



(l-z)" 

+ 



2 (l-z)l/2 2 (l-z)3/2 





2 

1 



2/(l-z)i' + ^(l-z)2 
1/2 
^(l-z)3/2 

^(l-z)»+l 

1 I 3 1 

4 (l_z)-V2 ^ 4-^ (l-z)5/2 
^ (a-l)a , „ a(a + l) 
^(l-z)<»+l (l-z)<'+2 



for m > 2 

(l-zl^ + l 
^(m-1)! 



^n-z) 

11/2]^ 



a; 7 



_ J^1/2U I ,,. [1/21„ 
1-1 -1a + m-l ' ^ 



' (l-z)''+"' 
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More detail discussions and related physical systems of these distributions will be given and discussed in the 
next section. Brief discussion about the weight Xk of each model follows. All the distributions, except the Poisson 
distribution (P), considered in Table I have several factors in the weight Xk- One factor in Xk is z'^ which is a 
k dependent geometric term and comes from assigning the same weight z to each constituents independent of the 
cluster or cycle classes it belongs to. Another factor of weight is independent of k such as the x in Table I which 
comes from assigning the same weight x to each cluster or the cycle class as a whole independent of its internal 
structure. These two factors, xz'^, are multiplied by a k dependent or independent prefactor. A geometric (Geo) 
distribution follows when there is no other weight factor beside x and z, i.e., no k dependent prefactor so that 
Xk = xz^ . The Geo with z = 1 for a finite N and with z = for fc > iV is the same as the Poisson distribution; 
both have /o = Nx. The negative binomial (NB) which appears frequently in various studies has a weight factor 
assigned to a cluster or cycle class given by xz'^/k. This has an extra size dependent factor of 1/k compared to the 
geometric distribution. The signal/noise model (SN) has a two part structure and interpolates between a Poisson 
and NB distribution. The geometric distribution is the signal component of SN while the NB distribution is the 
noise component of SN. The Lorentz/Gatalan model (LC) has in its weight a shifted Catalan number divided by 

22(/s-i)^ that is ^^^fcf"' = ^ ^ ^^k- ) ' f^'Ctor which is the weight for Geo model. The Catalan 

numbers given by ( ^ j /{k + 1) are 1, 2, 5, 14, • • • for = 1, 2, 3, 4, • • • and the shifted Catalan numbers given by 

/k are 1, 1, 2, 5, 14, • • •. The importance of this factor is shown in Fig. 2 of section HE. 

Of all these distributions, the NB has been the most frequently studied. Ref. [33] gives several sources for its 
origin. These sources include sequential processes, self-similar cascade models and connections with Cantor sets and 
fractal stucture, generalizations of the Planck distribution, solutions to stochastic differential equations. Becattini 
[34] have shown that the NB distribution arises from decaying resonances. The a. model of Ref. [8], which is a self 
similar random cascade process, leads to a NB like behavior. The stochastic aspects of the NB distribution have been 
discussed by R. Hwa [35]. Hegyi [36] has discussed the NB distribution in terms of combinants. As already mentioned, 
the LC model can be connected to a Ginzburg-Landau approach and also has an underlying splitting or branching 
dynamics and cascade like features. 

As can be seen from the arguments of the hypergeometric fimction F{a, 6; c; z) in Table I, the hypergeometric model 
with 6=1 and c = 2 (HGa) include Geo, NB, SN, LC as a special case of HGa depending on the value of a. Other 
models listed in Table I are based on random walks. The use of random walk results was originally due to Feynman 
[37] in his description of the phase transition in liquid helium. The random walk aspects arise when considering the 
closing of cycle of length k. We include them for completeness. Since the random walk in 1-dimension (RWID) is 
the same as LC except the missing 1/k dependence compared to LC, RWID can be extended to a generalized RWID 
(GRWID) similar to the generalization of LC to HGa. A random walk model in 2-dimension has an extra factor of a 
shifted Catalan number and k2^^^''^^^ factor compared to RWID and can also be generalized to GRW2D. 

Since fc = is excluded here, the partition function for these models are given simply by a hypergeometric function 
as Z = exp[^^^ Xfc] = exTp[xzF{a,b;c;z)] with various choices of a, b, c. For example the LC model has /o = 
^k^k = xzF{l/2,l;2; z) and the NB model has /o = xzF{l,l;2; z). The geometric model Xk = yz'' has /o = 
yzF{a, 1; a; z) = yzF(2, 1; 2; z) while the SN model is a combination of the geometric plus NB cases. These functions 
are special cases of /o = xzF(a, 1; 2; z) of the HGa model. The generalized random work in 1-dimension (GRWID) 
has /o = xzF{a,b;b; z) and the generaized RW in 2-dimension (GRW2D) has /o = xzF{a,a;l; z). The factorial 
cumulants for these models are summarized in Table II. Several cases with c = 2 have a canonical partition 
function Z„ which can be writen in terms of confluent hypergeometric functions U {u, v; w) and standard factor z"/n\ 
[32]. 



G. Generalized model of Hypergeometric (HGa) 

We consider in more detail the hypergeometric model with b — 1 and c = 2 (HGa) here since it includes the NB, 
Geo, and LC models as special cases. This generalized model is related with the hypergeometric function with 6=1 
and c = 2 with an arbitrary value of a, i.e., F{a, 1; 2; z), and has the weight of 

k Mfc-i z^ rja + k^l) 

"^ = ""^ = "¥ r(a) 

The asymptotic behavior of Xk is 
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Xk = xz^'k"-'^ /T{a) = xz'^k-'' IT{a) 
for large k using Stirling approximation. Thus 

T = 2-o 



(55) 



(56) 



which connects the parameter o to the physical Fisher critical exponent t for the HGa class of distributions. Its 
associated grand canonical partition function 



1 



is shown in Table I. Prom Table II, we have 



fm{x,z) = [a]. 



xz" 



(1 - ^)a+"»-l 



: X 



(a-1) \{l-zY''-^) 
r(a + m-l) ^™ 



[a m-l 



x" 



r(i-^) 



(o-l)(m-l) _ 



r(a) (1 - ^)«+™-i 
r(a + m- 1) /(l-^)(°-i)^ ""^ 



Via) 



a 



m-l = ^mK^ ^{x,z) 



Tia + m.-l)K^-\x,z) _ , 



(57) 



(58) 



(59) 



The normalized factorial cumulant, i.e., the reduced cumulant shows the hierarchical structure of HGa at the 
reduced cumulant level with Am = a~^'"~^-'r(a + to — 1)/T{a) where Km is related to K2- This property was realized 
for the NB distribution in Ref. [28] which is obtained for Eq.(59) with a = 1 giving = (m — 1)!. The result of 
Eq.(59) is a generalization of the NB result. 
Some moments for HGa are 



{A){x,z) = fi{z,x) 



xz 



(l-z)" 



X{x,z) = 



(l-z) 



fi{x,z) (1-a) z 



^(x,.) = K2(a.,.) = #^ = ^(l-.)(«-i) 



Since these relations give 



fl{x,z) X 



^{x, z){A){x, z) = K2{x, z){A){x, z) 



f2{x,z) 



fl{x,z) (1-2) 

the void parameters can be obtained in terms of the normalized fluctuation ^ and the mean number {A) = Ahy 



ar(l,0 



a^iA a + hi A /z + aA 



-(l-.)« = ^ 

z^ ' ^\a + iA 



(a-l) 



aA^ 



aA 



/2 \aA + f2 



a-l 



/o(i,0 = log^(AO = 



X 



/o 



a-l 
1 a 



a 



- 1 



X 



1 + 



l-a 



A (l-a) 

fm{A,o _ na+m- i) a 

A^ r(a) \a 



- 1 



a-l 

1 aA 
(T^X 



1 + 4 
aA 



a-l 



1 + 



aA 



- 1 



(60) 
(61) 
(62) 

(63) 

(64) 
(65) 
(66) 

(67) 
(68) 



for a given mean value of {A) — A and the fluctuation ^ or ^A or /2. Using these x(A,^) and z{A,^) we can easily 
find the multiparticle distribution Pa{A, ^) for a given values of mean A and fluctuation ^ using the recurence relation 
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of Eq.(6) for Za{A,^). We explicitly compare the distribution Pa for different values of a with given values of A and 
^ in Sect. III. 

Tabic I shows that the generalized HGa model becomes the Lorentz/Catalan (LC) model with a = 1/2, the negative 
binomial (NB) model with a = 1, and geometric (Geo) distribution with o = 2. However the NB should be considered 
as a a ^ 1 limit of HGa; 

lim /o = -a;log(l - z) (69) 

a— >1 

lim Z = {1- (70) 

a— >1 

lim/o(AO=^ln(l + ^^) =xln(l + /2/A) (71) 

a— >1 

limx(Ae) = 7^1og(l + e^) (72) 



Thus aNB distribution has x = ln(l +C^)/(CA) with x = l/^ while the LC distribution has x = (vl + 2^- 1)/(^A) 
with X = y/l + 2^A/{2^) and a Geo distribution has x = 2[1 - (1 + ^A/2)-'^]/{^A) with 1/x = ^(1 + iA/2)/2 as 
limiting expressions of a more general x given by Eq.(67) for HGa. Eq.(65) shows that the fluctuation ^ depends on 
the mean A for given x except for NB. Since > for ak = k with non-negative Xk, £,{A,x) = 1/x for a = 1 (NB), 
^{A,x)_= [A + + 4a;2]/(4a;2) for a = 1/2 (LC), ^{A,x) = [-x_+ v^x^ + 4:xA]/{xA) for a = 2 (Geo). However for 
small A and large x the fluctuation ^ is inependent of the mean A and ^{A, x) « a/x for HGa. 

The void parameters are ^ = and x = 1 for Poisson distribution. The reduced cumulants arc Km = for m, > 2 for 
Poisson distribution. The void parameters are C = ^ ^-iid x = ^"^^g^^"^"* for NB distribution. The m-th order reduced 
cumulant is Km = {tti — 1)!^™~^ for NB. For SN distribution, the void parameters are 



h ^ N {2S + N) 
f? x{S + Ny 
fo ^ xln{l + N/x)+S/{l + N/x) 
h S + N 



^~f^~x{S + N)^ ^ ' 

(74) 



with signal level S = ^^^^2 and noise lc;vel N = The SN model has important application to quantum optics 

and, in particluar, to photon counts from lasers [38]. Biyajima [39] has suggested using it for particle multiplicity 
distribution as does Ref. [33]. When the noise level N ^ 0, ^ ^ and x ^ 1- When the signal level 5 — > 0, ^ 1/x 

and X ^ f ln(l + N/x). The third order reduced cumulant is K3 = + ^) = ^# + N/xf{3S + N) for 

SN. For LC distribution, the void parameters are (, = and X = ^ \/l+'2^ - 1 with z = 2^^/(1 + 2^^). 

The m-th order reduced cumulant is Km = [■^]m-i{2£,)"''~^ for LC. The results of quantum optics, in the notation of 
Ref. [38], can be obtained [31] when x = Tn/2, z = 2W"//n'^ Q"^ = j'^ + 2W'y. The W is an integral of the Lorentzian 
line shape r(a;) = b/[{u! — ojq)^ + 7^]> T is the time, and 2a; Vl — z = ^T. 



H. Generalized random walk in one dimension 



The hypergeometric model with c — b has Xk = x z^ which is a generalized random walk process in one 
dimension (GRWID). For a = 1/2 the GRWID becomes the random walk in one dimension. Since F{a^b\c\ z) = 
F{b, a; c; z), it is easy to see that both HGa with a = 2 and GRWID with a = 1 are geometric models. For GRWID, 



IS 



we have /o = xjj^-^ which is /i of HGa as shown in Table II. The asymptotic behavior of Xk = Xj^^j^.z 

Xk = xz^k°'~^ /T{a) for large k and thus the Fisher critical exponent is r = 1 — a for the GRWID class of distributions. 
Its associated grand canonical partition function is 



^ [<Ak- 



z{x, z) = ef° = e^^'f'(«'^;^;^) = 



The factorial cumurants for GRWID are, from Table II, 



exp 



xz 



(1 - zY 



(75) 



/■ / \ _ [oi m I 

Jm{X,Z)— X^^_^^^^^_^Z -i-X 



[l-z) 
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\{a-l)z + m] 
Km{x, z) = X- r— — - [{a - l)z + m] 



(76) 



(1 - 2;)a+™ 
(1 - z)'' 



{1-z) 



o+l 



[{a-l)z + l]- 



x[{a- l)z + 1] 



[a]m-i [{a - l)z + m] 



[(a - l)z + 1] 



[a + 1] 



m-2' 



[{a — l)z + m] 
[(a - l)z + 2] 



[(a - + 
a[{a- l)z + 2] 



-I m-2 



1^2 — 



(77) 



Here the coefficient Am that appears in the hierarchical structure relation of Eq.(77) depends also on z which means 
that Ajn depends on the mean A and the normalized fluctuation ^. For HGa A^ is independent of A and ^. 
Some moments for GRWID are 



{A){x,z)=A = fi{x,z) 

fo{x,z) 



xz 



X{x,z) = 



(1 - z)«+i 
(1-^) 



[(a -1)0+1] 



h{x,z) [{a-l)z + l\ 



^/ ^ .X f2(x,z) a, [(a — l)z + 21 

^{x,z) = K2{x,z) = ^^^^-il-zr ' ^ 



Since these gives 



(,{x, z){A){x, z) = K2ix, z){A){x, z 



fi{x,z) X 

f2{x,z) 



[(a-l)0 + l]2 
az [(a-l)z + 2] 



h{x,z) (1-0) [(a -1)2 + 1] 
we have, for a given mean value {A) = A and the fluctuation ^ or /2, 



z{A,0 
x{A,0 



{a - 2)iA - 2a ± a^J{iAf + A^A/a + A 

2(a- l)(^^ + a) 
A (l-z)"+^ 

7[(o- i)z + i] 

A ( 2{S,A + a) 



2 + iA 



a£,A + 2o2 T a^J{S,AY +AS,A/a + ^ 



a+l 



z \a£,A ± a\/(CA)2 + A$,A/a + 4 



2(a-l)(eA + a) 



/o(AO=^ 



= a;2: 



2(a- l)(CA + a) 



X(A0 



(l-2;)« + 2a? + aV(^^)2 + + 4 

^ ^ / ^A + 2a+ V(^^)2+4^^/a + 4 \ 1 = 
A \ ^A±^{^Ay+AiA/a + A J {a-1) 

fmjA, _ [ (g- 1)Z + TO 



(a - 1)0 + 1 



(a - l)z + 2 



rn-2 



2 + ^A 

m—1 



(o + 2(m - 1))^A + 2(to - l)a ± 0^/(^-4)2 + 4^A/a + 4 



(a + 2)^A + 2a± a^y (^A)"^ + A^A/a + 4 

a^A±a^{^A)^ +4^A/a + 4: \ """^ 
(a + 2)^A + 2a ± aV($^)^ + 4^^/a + 4 i 



[a]m-i 



(78) 
(79) 
(80) 

(81) 
(82) 



(83) 
(84) 
(85) 



(86) 



for 



Here the condition < < 1 determines the right one of ± sign. Note here that Km{A,(,) = [aj^-i 

HGa. For the GRWID, the normalized factorial cumulant Km{A,^) has an extra dependence of the order m which 
also depends on the mean A and the fluctuation ^. 
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III. COMPARISON OF PROBABILITY DISTRIBUTIONS WITHIN HGA 



The hypergeometric case (HGa) which we studied in more detail in Sect. II G included various cases of Table I; 
Geo, NB, LC which arc distinguished by one parameter a of the hypergeometric model. Thus wc use HGa to compare 
various models for pion distribution and other particle distributions, with particle number n = /eri/j. 



A. Voids and void scaling relations 

Void analysis looks for scaling propoerties associated with x = .fo/(^); specifically, x is a function of the combination 
^(n) = j2/{n) where ^ = K2 is the coefficient of (n)^ in the fluctuation cr^ = (n^) — {n)^ = (n) For generalized 

HGa model % is given by 

^ (l + ^inVa)^-"-! 
^ {l-a)^{n)/a ^ > 

We show the void variable x as a function of ^(n) in Fig. 3. This shows that we can vary a to fit data. Ref. 
[28] claims the NB distribution (HGa with a = 1) fits reasonably well the void distribution for single jet events in 
e+e" annihilation but Fig. 3 shows all of the curves might fit such data up to ^(n) ~ 1 since the various curves are 
reasonably close up to this ^(n). Thus further investigations of this data is required to distinguish various models. 
Higher moments or cumulants might need to be compared for this purpose. Within HGa, K^{n, ^, a) = (^-^) for a 
given value of (n) = n and ^ according to Eq.(68). Thus A3 = m/^^ = (a + l)/a = 3, 2, 3/2, 4/3, 5/4 for a = 1/2, 1, 
2, 3, 4 independent of ^. The K3 may help in distinguishing various models for the data with a small value of ^, i.e., 
the region of < 1 in Fig. 3. The differences of Am between models with different value of a becomes larger as the 
order m of the reduced factorial cummulant becomes higher. 

Table HI shows the values of ^ and K3 of charged particle multiplicity distribution of jets for L3 data [40] and HI 
data [41]. From these values A3 = 4.1686 and —6.2653 for L3 and HI respectively. The corresponding values of HGa 
model parameter a are 0.31559 and —0.13764 respectively and the corresponding values of GRWID model parameter 
a are 0.14006 and —0.8516. The negative value of a makes the corresponding z value in HGa larger than 1 which is 
not allowed because of Eq.(65) for real x. Thus the HGa model cannot fit HI data with the same values of mean, 
fluctuation and K3 at the same time. The negative value of kz or ^3 causes the negative values of a, x and z in 
GRWID which are acceptable in Eq.(83) for real x. This model with a = —0.8516 fits HI data best for all the values 
of mean, fluctuation, K3, and x = /o//i- The Fisher exponent for this parameter is t = 1 — a = 1.85. The physical 
meaning of the negative x and z should be studied further similar as studied in Ref. [21]. Both HGa and GRWID 
models with small values of a (0.31559 and 0.14006 respectively) fit the L3 data well with the same mean, fluctuation, 
and the third cumulant at the same time. The corresponding values of Fisher exponents are r = 2 — a = 1.68 and 
T = 1 — a = 0.86 respectively. If we can extract the value of r from data then we may distinguish which model fits 
the data better. 
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TABLE III. Various cumulants of charged particle multiplicity distribution of jets for L3 data (all events of e+e~) [40] and 
HI data (pseudorapidity range 1 < ?7* < 5 with 80 GeV <W < 115 GeV of e'^p) [41] compared with models having the same 
mean and fluctuation. 



model 


a 


X 


z 


(n) 




/t3 


X 


L3 data 








20.463 


0.044238 


0.008158 




HGa 


0.5 (LC) 


18.948 


0.64419 


20.463 


0.044238 


0.005871 


0.74726 




1.0 (NB) 


22.605 


0.47514 


20.464 


0.044238 


0.003914 


0.71208 




2.0 (Geo) 


31.123 


0.31159 


20.463 


0.044238 


0.002936 


0.68841 




0.31559 


18.007 


0.74150 


20.463 


0.044238 


0.008158 


0.77636 


GRWID 


0.14006 


20.905 


0.66088 


20.463 


0.044239 


0.008158 


0.78558 


HI data 








7.7210 


0.069186 


-0.000764 


0.73242 


HGa 


0.5 (LC) 


10.394 


0.51653 


7.7210 


0.069186 


0.014360 


0.82028 




1.0 (NB) 


14.454 


0.34819 


7.7210 


0.069186 


0.009573 


0.80122 




2.0 (Geo) 


22.814 


0.21079 


7.7210 


0.069186 


0.007180 


0.78921 




3.0 


31.244 


0.15115 


7.7210 


0.069186 


0.006382 


0.78470 


GRWID 


-0.85160 


-4.4150 


-0.81168 


7.7211 


0.069185 


-0.000764 


0.72375 



17 



B. Probability distribution 



Once we know the generating function Z{x) or fo{x) we may study the probabiUty distribution P„ using Eq.(9) or 
the recurrence relation Eq.(6); 



The Pn for various cases of Table I are shown in Fig.4 and Fig. 5 with the same mean value (n) = n and the same 

fluctuation (, = K2 = hjn^ for fixed (n) = 10. Fig. 6 shows KNO plots of {n)Pn versus n/{n) for fixed (n) = 10 
and 20. Fig.4 shows that the various models considered here have almost the same distribution for small fluctuation 

= 0.01) and in this case they are very similar to a Poisson's distribution. For ^ = 0.05 the models are similar to 
each other except for larger n/n even though they are different from a Poisson distribution. For larger fluctuations 
such as ^ = 0.5, the models have very different forms even if they have the same mean value and fluctuation. Fig. 5 
shows that the probability distribution of these models differ in their form for fluctuations larger than ^ w 0.2. Ref. 
[33] shows the total charge distribution in hadronic collisions. They flt the data in the range of n = 6 ~ 13 using a 
SN distribution. The SN is a mixture of NB (« = 1) and geometric (a = 2) which are the cases we have shown in the 
figures. To determine which model fits the data best, we need to know the exact value of ^ and even higher moments 
of the data beside the mean value n. 

Fig. 7 shows the multiplicity distribution for L3 data and HI data discussed in Table III. This shows that all the 
models having the same mean and fluctuation also have very similar flts to the data for n/{n) smaller than 2. For 
n/{n) larger than 2, the NB (thin solid curve) or LC (thick dashed curve) fits best the L3 data. But for n/ (n) smaller 
than 1, the GRWID (thick solid curve) and HGa with a = 0.31559 (dash-dotted curve) flt the L3 data best and only 
these two cases give the correct value of K3 for L3. We can also see that the GRWID (thick solid curves) fits best 
the HI data for whole range of n/ {n). Only this GRWID model fits the value of K3 for HI data correctly. Thus we 
need to evaluate higher order moments of data (or at least up to third order, K3) beside the mean and fiuctuation to 
understand the distribution of data and its underlying mechanism. One difference between L3 and HI data is that 
the HI data has much larger probability for small n/{n) compared to the L3 data. This large probability at small 
n/ {n) causes the negative value for the third cumulant /ts for the HI data. 




(88) 




(89) 



(90) 
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n/<n> 

FIG. 4. P„ for fixed (n) = 10 for a = 1/2 (solid line), 1 (dash), 2 (dash-dot), 3 (dot), 4 (dash-dot-dot-dot) and for ^ = 0.01, 
0.05, and 0.5 in log scale on the left and linear scale on the right. For ^ = 0.01 all distributions become very close to Poisson 
(thin solid curve) already. (Po becomes large for large ^.) 
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I I I I I I I I I I I I I I II I I I I I I I I I I I I I I 




n/<n> 

FIG. 5. P„ for fixed (n) = 10 and for a = 1/2, 1, 2, 3, 4. The value ^ = /2/(n)^ are shown in each figure. The various 
choices of a for each curve are given in the figure caption of Figs. 3 and 4. 
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n/<n> 

FIG. 6. KNO plot of {n)Pn for fixed {n) = 10 and 20 for a = 1/2, 1, 2, 3, 4. 



1 








t \ 13 data 


'C^ HI data 


0.1 

c 
Q. 
A 
c 

^ 0.01 










Xi 

E 10-= 
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10-* 




\ 


10"* 


M., 





1 2 3 1 2 3 



Multiplicity, n/<n> 

FIG. 7. Charged particle multiplicity distribution in jets of L3 data [40] and HI data [41]. The crosses are the data and 
the curves are the fits with HGa and GRWID given in Table III. The thick solid curves are the GRWID fit, the thick dashed 
curves are LC model, the thin solid curves for NB model, the thin dashed curves for Geo model. The dash-dotted curves are 
for HGa model with a = 0.31559 for L3 data and with a = 3.0 for HI data. 
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The KNO behavior for n ^ oo can also be studied for the models listed in Table I. According to KNO scaling, for 
distributions with a large mean value, the distribution (n)P„ becomes model independent in the new variable n/n, 
i.e., variable scaled by mean value. In general, any distribution becomes Gaussian for large mean in) according to 
the central limit theorem; specifically 



P„(n, a) 



1 



: exp 



(n — n)^ 



= 1 + 7^ 

in 



-1/2 



2a2 
1 



: exp 



1 /n \2 / 1 



in 



(91) 



with the mean 

^2 _ /^2\ _ /^\2 



= n and the standard deviation a which is related to the reduced factorial cumulant ^ by 
a" = {n") — [n)" = n + ^n^. This means that the KNO scaling follows when the fluctuation is given by cr^ = ^n^ with 
a constant ^ or when ^ 1. Thus to compare KNO scaling properties of nP„ for different models, the fluctuation 
of these models should have the same value of ^ = /2/ri^. For a small ^, the Poisson component of the fluctuation 
0-2 = n + becomes dominant and thus the KNO scaling behavior is broken. For large ^ the Poisson component is 
negligible and KNO scaling is realized, i.e., 



nPn{n,i) = 



: exp 



2^ \n J 



(92) 



The KNO plot of Fig. 6 shows that the various distributions have no KNO scaling property for small fluctuation 
{£, = 0.1) but show a KNO scaling property for large fluctuation (^ = 0.5). The effect of mean on nP„ is larger than 
the difference between different model for ^ = 0.1 while the effect of mean on nP„ is much smaller than the difference 
between different models exhibiting KNO scaling behavior for ^ = 0.5. Ref. [33] shows the total charge distribution 
in hadronic collisions studied using SN model. From their fits we can extract the corresponding fiuctuation which are 
^ = 0.05 ~ 0.5 and n = 6 ~ 13. This means that the KNO scaling behavior is marginal for these data, i.e., just fitting 
the distribution nP„ of the data does not show clear evidence for KNO scaling. We need to evaluate the explicit 
values of the mean number n and the fluctuation ^ to check the KNO behavior of this data; the value of should 
be large enough to show KNO scaling behavior as can be seen from Eqs.(91) and (92). Table III shows that ^ < 0.1 
and ^(n) < 1 for both data indicating that the KNO scaling behavior would not be clearly exhibited in jets for L3 
and HI. 



C. Sequential procedures and compound Poisson distributions 

A Poisson distribution plays a very important role in physics. As already noted, in statistical physics, Maxwell- 
Boltzmann statistics leads to Poisson probabilities. Other distributions are compared to the Poisson distribution 
which acts as a benchmark for comparison. The distributions considered in this paper can have large non-Poissonian 
fluctuations. The purpose of this section is to show how they can be rewritten as a compound process or sequential 
process involving one aspect that has a Poisson character. As an example the negative binomial distribution can be 
obtained from a compound Poisson- logarithmic distribution as discussed in Ref. [33]. Here, wc extended this result to 
include the other distributions and wc also show that the final distribution can be obtained from compounding it with 
another distribution, such as the negative binomial. In general, the underlying picture for a sequential process involves 
a two step procedure in which the observed particles arise from the production of "clusters" with the subsequent decay 
of each cluster producing its distribution of particles. The final distribution is obtained by compounding the probabiUty 
distribution of the clusters with another distribution coming from each cluster and suming over clusters. Specifically, 
the observed particles or members in system arise from production of M = c clusters with probability distribution Pg. 
This is sequentially followed by each cluster decaying into ka particles with the probability Pfe^ with a = 1, 2, • • •, c. 
The probability of observing n = J2k ^nu = particles is then obtained by a compound probability expression 

c 

= E E n (93) 

c {fc„} a=l 

A negative binomial distribution can be obtained when Pc = (c)'^e~^°^/c! and P^^ = (g'^^/fca)/ ln(l/p). Here c = 
M = Y^^Uk and (c) = (M), p = {I - z) = (1 + {n)/x)-^ and g = 1 - p =" z = {{n)/x)/{l + {n)/x). Also 
Nc = (c) = a;ln(l + {n)/x) and ric = {n)/Nc = ((n)/x)/ln(l + {n)/x). This structure can be generalized as follow. 
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Since the generator of Poisson distribution is an exponential, i.e., the expansion of exponetial gives the Poisson's 
distribution P^{n), 

^ = E ^ = ^^E^-^^ = E P'mW (94) 

M=0 ■ M ■ M=0 

The grand partition function or the generating function Z = eJ" for any distribution can be represented as a Poisson 
distribution whose mean value is the void variable /o or the grand potential fl = — /o- On the other hand and in 
general we can rewrite the void variable as 

/o(f)=lnZ(f)=^Xfe=AA^Pfe(f) (95) 

k k 

where M = fo = Y^^. Xk and 

VkiS) = ^ (96) 
The Vk{x) can be connected to its generating function Q{x,u): 

g{x,u) = ^(1 - ufVk{x) = ^Y.'^kil - uf (97) 

k k 



Thus the generating function G{x, u) of P„ can be expanded in terms of Vk as 

G{x, u) = ^(1 - uTPnix) = ^ E - = 



gA/'a(x,u) 



Z{x)^ "V - ' 

-, M 



.AT 



M=0 ■ M 



= E E IlKl - ")'^'^]"'' (98) 

M UM k 

where P^{Af) = e~^Af^ /Ml is the Poisson distribution with the mean of A/" = /o- Here the sum over um is the sum 
over partitions n with a fixed M = J2k '^fc- Thus we have 

P„(x) =^pP(AA)^n^r(^) (99) 

M UM k 

with M = Hk and n = /cnfc. Any distribution obtained from a generating function of the form of Z{x) = 
eJ° = cxp[^^ Xfe] can therefore be decomposed as a compound Poisson's distribution with some other distribution 
Pk = Xk/ fo obtained from the weight Xk- 

The sequential nature of a process is explicitly shown on Eq.(99). The observed particle multiplicity distribution 

arises from a two step process in which M = Uk clusters are first distributed according to a Poisson distribution. 
This is then sequentially followed by breaking each of the clusters of type k into k particles with probability 
'Pk = Xk/N and with n = X^fc^^^fc- The probability associated with a given M and n with x is PM{x,n) = 

pimYik-PT^pimiik^^k/Mr-. 

As an illustration we consider the LC model with Xk = xCkz'^ /2'^^'' ^\ Using the evolutionary variables [31] 
X = (3/Ap and z = 4p(l -p) then TV = ^^fc = /3 for p < 1/2 as already noted in Sect. HE so that Pj^ = e'l^fi^ /M\. 
The Xk = pCkp''~^{l — p)^ so that Tk = Xk/Af = Ckp''~^{l — p)'^- The underlying diagram associated with Pk 
are shown in Fig. 2. For a negative binomial (NB) distribution, Xk = xz^ /k and thus Pk = Xk/M is generated from 
N = "Ylik -^fc ~ ~^ \T^{\ — z). Therefore the NB is a compound Poisson-Logarithmic distribution as shown in Table IV. 

As another example, we consider the HGa model with general a instead of a = 1/2 for LC or a = 1 for NB. The 
weight Xk has the structure of the probability Pk{x, z) of NB distribution given by Eq.(90), i.e., 

2;''r(a + fc-l) X 1 „NR/ 

= 'kl Via) = ^ (13^^^ - 
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thus Pfe = [l-(l- 



LpNB HGa. Therefore the HGa P„ distribution is a compound Poisson-NB distribution. 



This may interpreted as a sequential process in which clusters with a Poisson cluster distribution Pc breakup into 
particles with a particle distribution Vk given by a NB distribution with parameter a. For the various models 
considered here with their Xk listed in Table I, the corresponding distribution Vk and the normalization factor A/" = /o 
are listed in Table IV. We can further see that HGa can be looked as a compound Poisson-Poisson-Logarithmic 
distribution, i.e., a distribution having three sequential steps; Poissonian breakup into clusters — > Poissonian breakup 
of each cluster logarithmic breakup of each of them. 



TABLE IV. Poissonian sequential distribution for various models of Table I. 



Model 



Weight N 



Distribution Vk 



Comments on Vk 



P 
Geo 

NB 
LC 

HGa 
GRWID 



(n) 

2x [i - yr^] 



(1 - z)z^-'' 

1/2 



(1-^)1 



zfc r(fc-i/2) 

k\ r(l/2) 



-az z_ 
JsL. 



Monomer only or Uniform for N species 
Uniform with constituents 
logarithmic with constituents 

NB with constituents with Vo ~ 

NB with constituents without A; = 
NB with X = a 
Poisson (exponential) 
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Because of a unique role played by the Poisson distribution and the form e^" = exp[^j, Xk] of the generating 
function, the cluster distribution P^, is usually taken to be a Poisson. However, as noted before, other divisions are 
possible. Using the same approach used above for Pc = P^, we can expand any distribution using a NB instead of 
Poisson, i.e., with = P^ using the form of the generating function Z{x,z) = (1 — ^)~^ for the NB. Replacing z 
by M{z), the normalization factor of a new distribution, we have 



The g{u) is 



while the G{u) is 



Z{x, z) = [l- AA(2)]-" = [1 - M{z)g{u = 0)]-" 
N{[1 - u]z) 



g{u) 



N{z) 



G{U) = ^(1 - uTPn = 



Z{x,[l-u]z) _ [1-Afg{u) 



Z{x,z) 



(101) 



(102) 



(103) 



Thus we have 



G{u) = £ i^-^r^^:^mur = £ pT{^,^fmu)r 



M=0 



M=0 



NB 
M 



M=0 



n M 



M=0 {nk}M k 



pux) = ^p^^{x,M) J2 n^r(^) 



(104) 
(105) 



M 



The result of Eq.(105) shows that the distribution P„ can be written as a compound probability distribution of a 
negative binomial P^ with another probability Vk distribution generated from g{u). For the case of Afiz) = e^, 
which may be considered as the fugacity = for a particle with chemical potential n = z, g can further be 
decomposed as 



g{u) = 



Af{{l-u)z) _ e^^-")^ 
Af{z) ~ 



= £(i-«)'=pp(^) 



(106) 



fc=0 



i.e., Vk for this case is Poisson P^{z). li M = 1 - p/"/^, then Vk = PlUo/x) without fc = and Z{x,z) = 
[1 - (1 - g/n/^)]"' = eJ«. For /o = Y.k 

Xk given in Table I, the Vk becomes the same probability with Xk replaced by 



can be decomposed as a sequential process 



Xk/x. As an example the HGa with Z{x, z) = exp ^y3^[1 — (1 — z)^ 

consisting of a NB distribution of clusters with Z{x,Af = 1 — e^°^^) which is then followed by a breakup of clusters 
distributed with a HGa distribution with Z{x = l,z) but without voids, i.e., with Vq = 0. Thus this decomposition 
separates the parameter x assigned to cluster from other parameters. 



IV. CONCLUSION 



Event-by-event studies from high energy collisions are being used to study the details of particle multiplicity 
distributions as, for example, those associated with pions. Such studies not only give information about the mean 
number of particles produced, but also information about fluctuations and higher order moments of the probability 
distribution which are important tools for studying the underlying processes and mechanisms that operate. They are 
also useful in distinguishing various phenomenological models. Issues associated with fluctuations play an important 
role in many areas of physics and departures from Poisson statistics are of current interest. One purpose of this 
paper was an investigation of various models of particle multiplicity distributions that can be used in event-by-event 
analysis. These various phenomenological models are developed using a general form of a unified model which is based 
on a grand canonical partition function and an underlying weight arising from Feynman's path integral approach to 
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statistical processes. A resulting distribution has three control parameters called a, x, z. The relationships of these 
parameters to various physical quantities are discussed. One important result is the connection of the parameter a 
to the Fisher exponent r; namely t = 2 — o for HGa. This connection arises from a parallel we developed between 
the model for particle multiplicity distributions considered here and our previous approach to cluster yields. Since an 
exact description of particle multiplicity distributions is not known, we have considered several cases with different 
r's or a's which are contained in our unified description. Moreover, many of the existing distributions currently used 
in particle phenomenology are shown to be special choices of r or a which appear in a generalized hypergeometric 
model called HGa. These include the Poisson distribution coming from coherent emission, chaotic emission producing 
a negative binomial distribution, combinations of coherent and chaotic processes leading to signal/noise distributions 
and field emission from Lorentzian line shapes producing the Lorentz/Catalan distribution. Using the HGa model 
combinants, cumulants and moments are discussed and a physical significance is given to combinants in terms of 
the underlying partition weights of a Feynman path integral approach to statistical processes. The parameter a or 
Fisher exponent r is shown to play an important role in the behavior of the combinants which manifest itself in 
various physical relationhips. The HGa model and its associated special cases are used to explore a wide variety 
of phenomena. These include: linked pair approximations leading to hierarchical scaling relations on the reduced 
cumulant level, generalized void scaling relations, clan variable descriptions and their connections with stochastic 
variables and branching processes, KNO scaling behavior, enhanced non-Poissonian fluctuations. Models based on an 
underlying random walk description are also discussed. 

In this paper we compared various particle multiplicity distributions within the hypergeometric model HGa. Our 
results show that even though various distributions have the same mean and fluctuation, the distribution itself or 
the underlying mechanism could be different. Comparisons within the HGa model also show that just comparing 
void variables x and ^ or mean n and fluctuation /2 or a is not enough to distinguish different models that descibe 
particle multiplicity data. Thus, to find the correct distribution and underlying mechanism from various data more 
information than just the mean and fluctuation are necessary and new variables should be found which are quite 
different between different models. For example, higher order reduced factorial cumulants need also to be evaluated 
such as the third order cumulant K3. 

Applications of our approach to the charged multiplicity data of L3 and HI are given. The mean (n) , fluctuation 
variable ^, void scaling variable Xj and third order reduced cumulant variable K3 obtained from these experimnets are 
compared with various models discussed in this paper. The value of the third cumulant Ks shows that the GRWID 
model can fit the charge multiplicity data of L3 and HI. The HGa model also gives a fit to the L3 data, but is not 
as successful as GRWID with respect to the third order reduced cumulant of the HI data. Both have very similar 
multiplicity distributions as shown in Fig. 7 which compares the two models with the charged particle multiplicity 
distributions of jets for L3 and HI data. Differences appear in the high multiplicity events with reduced multiplicity 
n/(n) > 2.5. The reduced third cumulant is negative for the HI data of Table HI. Negative values of a cumulant 
require negative values of the parameter a which are allowed in GRWfD if the parameters x and z are also taken 
as negative, but a negative value of a is not acceptable in HGa since this would require a value of z greater than 1 
and thus a complex x. We also discuss the KNO scaling behavior for these data and conclude that the KNO scaling 
behavior would not be clearly exhibited in the jets for L3 and HI. 

In this paper we have also generalized the compound distribution that arises from sequential process which may 
reveal the dynamical structure of the distribution. Specifically, the underlying sequential picture involves a two step 
process where the final distribution arises from the production of clusters followed by a subsequent decay of the 
clusters. For the HGa model, the final distribution is obtained from compounding a Poisson distribution of clusters 
with a NB distribution coming from the decay of each of the clusters. The HGa may arise through a three step 
sequential process of Poisson-Poisson-Logarithmic compound distribution. It is also shown that the HGa can arise 
from a two step sequential process of a NB distribution followed by a new HGa with a different mean value. 

This work was supported in part by Grant No. 2001-1-11100-005-3 from the Basic Research Program of the Korea 
Science and Engineering Foundation and in part by the DOE Grant No. DE-FG02-96ER-40987. 
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